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ABSTRACT 

The hard X-ray lightcurves exhibit delays of ~ 1 s with respect to the soft X-ray 
lightcurves when the microquasar GRS 1915+105 is in the state of frequent, regular 
outbursts (states p and k of Belloni et al. 2000). Such outbursts are supposed to be 
driven by the radiation pressure instability of the inner disc parts. The hard X-ray 
delays are then caused by the time needed for the adjustment of the corona to changing 
conditions in the underlying disc. We support this claim by the computation of the 
time evolution of the disc, including a non-stationary evaporation of the disc and mass 
exchange with the corona. 
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1 INTRODUCTION 

The observations of Galactic black hole binaries (GBH) im- 
ply the coexistence of a relatively cold, optically thick ac- 
cretion disc, responsible for a thermal disc-blackbody com- 
ponent in their soft X-ray spectra, with a hot, optically thin 
medium that is the source of power-law spectral tail in the 
hard X-ray band (see e.g. Done 2002 for a review). At least 
in the High and Very High spectral states, the latter may 
have the form of a corona above the accretion disc, which 
means that the geometrical configuration of these two me- 
dia is such that both are extending down to the last sta- 
ble orbit around the black hole, being vertically separated 
from each other. The hard X-rays are produced via Comp- 
ton upscattering of seed photons coming from the under- 
lying disc, which basically requires the energy dissipation 
within the corona. Apart from this radiative coupling, also 
a mass transfer, i.e. evaporation and condensation of matter 
between the disc and corona is possible (Meyer & Meyer- 
Hofmeister 1994; Rozanska & Czerny 2000). 

Radiation pressure instability (e.g. Taam & Lin 1984; 
Lasota & Pelat 1991) seems to be a plausible mechanism 
to account for the characteristic variability of the Galactic 
microquasar GRS 1915+105 (Mirabel & Rodriguez 1994; 
Taam, Chen & Swank 1997; Belloni et al. 2000). Exemplary 
lightcurves were analyzed recently by Naik et al. (2002). 
Some of these lightcurves exhibit a very regular shape of 
outbursts, that can be possibly related to the disc variabil- 
ity, while others are more chaotic and probably driven by 
other mechanism, e.g. jet emission. Time-dependent accre- 
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tion disc model with jet emission was studied e.g. by Nayak- 
shin et al. (2000) and Janiuk et al. (2002) . Recently, Watarai 
& Mineshige (2003) analyzed the oscillations of this source 
allowing for the temporary evacuation of the inner disc. 

If the disc instability is the primary cause of regular 
outbursts, we may expect that the hard X-ray coronal emis- 
sion would be significantly delayed with respect to the soft 
X-ray disc emission. In the present paper we check the hy- 
pothesis of the radiation pressure instability as the outburst 
driver by analyzing the observed time delays between the 
soft and hard X-ray emission. We study the behaviour of 
GRS 1915+105 in various variability classes. We compare 
the observed time delays with the exemplary results of our 
theoretical model. In this model we compute numerically 
the time evolution of the disc-corona system with the mass 
exchange between the two media. 

The structure of this article is as follows. In Section 
we present the results of analysis of the microquasar GRS 
1915+106 observations, obtained by the Rossi X-ray Timing 
Explorer (RXTE). The X-ray lightcurves were examined to 
determine the time lags between hard and soft X-ray bands 
in different modes characteristic for the variability of this 
source. In Section|2]we describe our model and assumptions 
about the disc and corona structure. The initial state of 
the disc plus corona system, as well as the prescription for 
the mass exchange, are described in Section 13.11 while the 
time-dependent equations, according to which this system 
subsequently evolves with time, are given in Section ^. 21 The 
results of the evolution are given in Section |1] We discuss 
our model and compare its predictions with observations in 
Section The conclusions are given in Section |S| 
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Figure 1. The X-ray lightcurves of GRS 1915+105 obtained from 
RXTE observation from June 20 2000 (class p), in the energy 
bands 1.5-6 keV and 6.4 - 14.6 keV. 

2 TIME DELAYS 

Fourier resolved time delays in the lightcurves of GRS 
1915+105 were analysed by Muno et al. (2001) for the 
periods of extended hard steady state. In the present pa- 
per we concentrate mostly on periods with significant out- 
bursts. In this case Fourier resolved phase lag approach is 
not best suited since clear semi-periodic signal dominates 
each lightcurve. Therefore, we restore to the simplest direct 
delays, as measured by the cross-correlation function. 

In order to determine the time lag between soft and 
hard X-ray emission we performed the Fourier analysis of 
these lightcurves by means of the Fast Fourier Transform 
(FFT) method. The cross-correlation function of two peri- 
odic functions F(t) and G(t) is defined as: 

Corr(At) = J F(t)G(t + At)dt. (1) 

For our analysis we select lightcurves representative for 
various variability classes, as studied in detail in Belloni et 
al. (2000). We choose the exemplary observations of GRS 
1915+105 made between 1996 and 2000, available through 
the public RXTE archive. The data were binned to 0.256 
seconds and the lightcurves were generated for the two en- 
ergy bands separately: 1.5-6 keV (PHA channels 0-14) and 
6.4-14.6 keV (PHA channels 15-35). Each lightcurve consists 
of one or more intervals, and in the subsequent analysis we 
compare the single intervals between each other. 

In Figure we show an exemplary lightcurve of the 
microquasar observed by RXTE on 20 June, 2000. This 
variability pattern belongs to the class p of Belloni et al. 
(2000). The cross-correlation function is calculated for the 
two lightcurves obtained in this observation and shown in 
Figure|2| The maximum of this cross-correlation function de- 
fines the time lag, At, between the two lightcurves. In this 
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Figure 2. The normalized cross-correlation function of the two 
lightcurves, 1.5-6 keV and 6.4 - 14.6 keV, shown in Fig. The 
main peak is shifted to -0.768 s, which means that the hard X-ray 
lightcurve lags the soft X-ray one by 0.768 s. 



Table 1. Time lag of the hard X-ray lightcurve (6.6-14.6 keV) 
with respect to the soft X-ray lightcurve (1.5-6 keV). Ati is cal- 
culated for the non-smoothened lightcurves, binned by 0.256 sec, 
and At2 is calculated for the lightcurves smoothened over AT = 1 
s. 



PID 


Date obs. 


Class 


At! [s] 


At 2 [s] 


20402-01-28-00 


18/05/97 


a 


0.0 


0.0 


20402-01-41-00 


19/08/97 


S 


0.0 


0.0 


20402-01-37-00 


17/07/97 


7 


0.0 


0.0 


20402-01-44-00 


31/08/97 




0.0 


0.0 


10408-01-15-00 


16/06/96 





0.0 


0.0 


20402-01-36-00 


10/07/97 


A 


0.0 


0.256 


20402-01-33-00 


18/06/97 


K 


0.256 


0.512 


50703-01-15-01 


20/06/00 


p 


0.768 


1.024 



case it is equal to At — 0.768 s. This means that the hard 
X-ray lightcurve lags the soft X-ray one by 0.768 seconds, 
which in this case is roughly 1.5% of the outburst duration. 
Other results are given in Tabled 

Regular, large timescale outbursts are characteristic for 
the lightcurves of type p and k. In other cases, the rapid, 
stochastic variability is overimposed on the more regular 
outburst pattern of a longer timescale (classes A, /3, and 9) 
or even dominates the whole observation interval (classes a, 
5 and 7). This rapid variability does not exhibit any time 
lags between hard and soft X-ray bands and is probably not 
correlated with the disc/corona system instability. There- 
fore, in order to determine more accurately the correlation 
between hard and soft X-rays in case of long outbursts, we 
filter out the rapid, small scale variations. We smoothen the 
lightcurves by applying the running mean filter (Stull 1988), 
to remove the variations shorter than (arbitrary) AT = 1 
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sec. The resulting time lags are given in the last column of 
Tabled We see that the time lags between regular outbursts 
are even more pronounced, if these outbursts are cleared 
from the stochastic variability. On the other hand, the time 
lag Af2 remains only marginally detectable, or equal to zero, 
in the lightcurves that exhibit mostly the stochastic varia- 
tions. 

The measured lags calculated for the two latter obser- 
vations, k and p, seem to be exceptionally high for an X- 
ray binary (Zycki, private communication). However, this is 
not surprising since only GRS 1915+105 exhibits outbursts 
while other sources show only stochastic type of variability, 
which must be of different nature. 



3 TIME-DEPENDENT DISC/CORONA 
MODEL 

Accretion discs are known to be locally unstable in cer- 
tain temperature (corresponding to the accretion rate) and 
density ranges, in which the cooling and heating balance is 
strongly influenced either by atomic opacities or by the ra- 
diation pressure. Such instabilities do not disrupt the disc 
completely, but lead to repetitive outbursts. 

The first type of instability, connected with the partial 
hydrogen ionization (Smak 1984; Meyer& Meyer-Hofmeister 
1984), operates in accretion discs of binary systems in the 
range of radii of the order of 10 4 — 10 5 i?schw and is responsi- 
ble for the luminosity changes on timescales of months, e.g. 
Dwarf Nova outbursts. 

The other instability, caused by the radiation pressure 
domination over the gas pressure, operates in the innermost 
regions of the disc, where the temperature exceeds ~ 10 6 
K, and is responsible for the disc variability on the shortest 
timescales (of the order of tens - hundreds of seconds). This 
radiation pressure instability was first noticed in Pringle, 
Rees & Pacholczyk (1974) and studied in Lightman & Eard- 
ley (1974) and Shakura & Sunyaev (1976). Abramowicz et 
al. (1988) found that radial advection has stabilizing effect 
on the disc at high accretion rates and the time-dependent 
calculations of the disc limit-cycle behaviour were performed 
by Honma et al. (1991) and by Szuszkiewicz & Miller (1998). 

The innermost disc regions in which the radiation pres- 
sure instability is possible, are covered by the hot corona. 
Direct comparison of the extension of the radiation pres- 
sure domination zone and the corona covered zone, resulting 
from the model of the stationary, two-temperature corona, 
was performed by Janiuk & Czerny (2000). In the time- 
dependent model of radiation-pressure instability proposed 
in Janiuk, Czerny & Siemiginowska (2002) we used a simpli- 
fied description of a stationary corona above the fluctuating 
disc, parameterized by a constant value of the fraction of 
gravitational energy dissipated in the corona (/cor)- Here we 
extend our model to the case of a non-stationary corona, 
that forms due to the continuous evaporation of material 
form the disc surface. Due to the mass exchange with the 
disc the corona follows its time-dependent behaviour and 
therefore periodic changes are expected also in the hard X- 
ray luminosity. 



3.1 Assumptions and model parameters 

3.1.1 Disc 

First we describe the initial steady state disc model, from 
which we start our subsequent calculations. Throughout the 
calculations we use the vertically integrated equations of the 
disc structure, as the disc geometrical thickness H is always 
small (H/r ~ 0.01 in the quiescent disc and H/r ~ 0.1 in 
the outburst; see Section [4.11 . 

The ang ular velo city of the disc is assumed to b e Ke p- 
lerian, fi = yj GM/r 3 , and the sound speed is c s = y/P/ p = 
QH. Here p is the gas density in g/cm 3 , M is the mass of the 
accreting black hole and G is the gravitational constant. A 
non-rotating, Schwarzschild black hole is assumed and the 
inner radius of the disc is always at 3 Jischw The outer ra- 
dius is equal to 300 -Rschw, and at this radius a constant mass 
inflow, parameterized by the external accretion rate M BX t is 
assumed. Only the innermost zone up to ~ 100i?schw is the 
subject to radiation pressure instability, while the rest of the 
disc is stable (the exact value of the radial extension of the 
unstable zone depends on the model; see Section f4.lt . The 
mass of the black hole is assumed to be IOJWq. 

For the disc heating we assume that the viscous stress 
tensor is proportional to the total pressure P: 



-aP, 



(2) 



and the vertically integrated heating rate is 

Q+ sc = ^aSlHP 

where a is the viscosity parameter given by Equation!?] The 
total pressure P consists of the gas and radiation pressure: 



(3) 



+ P* 



P e ,s = pT 



(4) 
(5) 

(6) 



where T is the mid-plane temperature, and k, m p and a are 
physical constants. 

The angular momentum transport in the accretion 
disc is driven mostly by the magneto-rotational turbulences 
(see Janiuk et al. 2004 and references therein). Since the 
magnetic fields are quickly expelled from the radiation- 
dominated disc (Sakimoto & Coroniti 1989), the transport 
efficiency, and in turn the viscosity parameter, must de- 
crease when the radiation pressure becomes dominant. On 
the other hand, when the contribution from radiation pres- 
sure is only moderate, it may still be possible to couple the 
radiation to the particles. Here we adopt a modified viscosity 
law for the accretion disc (Chen & Taam 1993; Nayakshin, 
Rappaport & Melia 2000): 

(1+C/Co) 



a = a.® 



(i + mo) 



(7) 



where £ = -Prad/fgas (for the discussion of this parameter- 
ization see Section 15.111 . For the model parameters ao and 
£o we assumed the values 0.01 and 8.0, respectively. This 
prescription implies that for small to moderate values of £ 
we have effectively the disc heating proportional to the total 
pressure, while for large values of £ the viscosity is propor- 
tional to the gas pressure. Therefore the radiation pressure 
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instability may still operate, contrary to the so-called /3-disc 
prescription (Lightman & Eardley 1974). 

The cooling in the disc is due to advection and radiation 
and the radiative cooling is equal to: 



(8) 



where r is the optical depth, E = pH is the gas column 
density in g cm -2 , c and a are physical constants, and we 
adopt the electron scattering opacity k = 0.34 cm 2 /g. 

The advective cooling in a stationary disc is determined 
from the global ratio of the total advected flux to the to- 
tal viscously generated flux (e.g. Paczyriski & Bisnovatyi- 
Kogan 1981; Muchotrzeb & Paczynski 1982; Abramowicz et 
al. 1988) 



Qadv — 



and 



2rPg adv 
" 3pGM 



g adv = (12 - 10.5/3) 



dlnT 
d In r 



(4 - 3(3) 



d\np 
d In r 



(9) 



(10) 



Here j3 is the ratio of the gas pressure to the total pressure 
P = Pgas/P = 1/(1 + £). In the initial stationary disc we 
assume that q a( jv is approximately constant and of the or- 
der of unity (in the subsequent evolution the advection will 
be calculated more carefully, with appropriate radial deriva- 
tives) . 

In order to calculate the initial steady-state configura- 
tion, we solve the energy balance: F to t = Q^ isc = Q adv + 
Qrnd- H ere the total energy flux ftot dissipated within the 
disc at a radius r is calculated as: 



3GMM .. . 
Ftot = — — —f(r) 



87T7" 3 



where 



/(r) = ^ _ ( 3^chw )3/2 \r 



Rsc 



2R Sc 



(11) 



(12) 



is the boundary condition in the pseudo-Newtonian poten- 
tial (Paczynski & Wiita 1980). We choose the initial accre- 
tion rate M that is constant throughout the disc and is low 
enough for the disc to be on the stable gas-pressure dom- 
inated branch, where neither radiation pressure nor advec- 
tion is important (see Figure [3J . The initial model is cal- 
culated by means of a simple Newtonian method, through 
which we determine the radial profiles of density and tem- 
perature, as well as the disc thickness. 



3.1.2 Corona 

The corona is assumed to be geometrically thick and opti- 
cally thin. Therefore its height is assumed to be equal to 
radius, H cov — r, and the pressure in the corona is only gas 
pressure due to ions: 



Pcot — Fgas — PcorTct 



(13) 



The contribution from the electrons is neglected. The corona 
is hot and its ion temperature is assumed to be equal to the 
virial temperature: 



GM m D 



(14) 



The initial configuration of the corona is computed un- 
der the assumption that its optical depth is equal to unity: 



— ftS cor — 1.0 



(15) 



and therefore the corona has a uniform surface density. 
When the time evolution starts, the proper solution devel- 
ops in the middle parts of the disc, but will be fixed at the 
boundaries R ln and i? ut by the above condition. 



3.1.3 Mass exchange (prescription I) 

The mass exchange rate in the vertical direction, between 
the disc and corona, is equal to the ratio between the locally 
generated flux used to evaporate the disc material and the 
energy change per particle: 



Fm v 



AE/m p fcT c( 



(16) 



(measured per surface unit, g/s/cm 2 ). 

In a stationary disc the generated flux depends on the 
accretion rate and the disc radius (see Eq. II II . We as- 
sume that during the time evolution the energy flux lead- 
ing to evaporation preserves this dependence. Since both 
the energy dissipated within the corona and within the 
disc can lead to disc evaporation, we assume that the en- 
ergy flux is proportional to the sum of the disc and corona 
accretion rates, taken with different numerical coefficients: 
F oc (0.5BiM C or + -B2Mdi S c). The coefficients Si and B 2 are 
in the range from to 1 and express the fraction of the gen- 
erated flux that is used to drive the evaporation. The share 
of the corona is always lower than the half of the total corona 
flux, since half of the flux from the corona is directed toward 
the observer, whereas the other half is directed towards the 
disc and there reprocessed. 

The total accretion rate is M = M cor + Mdisc and may 
be locally constant in case of a stable disc. When the disc 
is unstable, both local accretion rates strongly depend on 
time and radius, and the relative contribution of the disc 
and corona to the total flow also vary. 

Expressing the accretion rates through the local vari- 



ables Mcor = 27rE 

c 



and A/disc = 27r£dis 



c rvf lsc in the 



corona and disc respectively, we obtain a useful formula for 
the vertical mass transfer: 



4r 



/(r)(0.5BiE c 



■ t> v disc\ 



(17) 



This formula is qualitatively similar to m z ~ Afc/'/r 3 ' 2 , de- 
rived by Rozahska & Czerny (2000) under specific assump- 
tions about the disc/corona coupling. However, their model 
described a stationary evaporation and did not depend on 
time. In our case the corona has to develop above the disc 
from the initial, uniform, very low density distribution, and 
this rise is provided by the sum of disc and corona accretion 
rates. (In the above, we denoted the quantities in the disc 
with the subscript "disc" , in order to discriminate them from 
the coronal ones. Please note, that whenever the physical 
quantities appear without any subscript, they also refer to 
the disc, and the coronal quantities are always distinguished 
with the subscript "cor".) 
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3.2 Time evolution 

Having computed the initial disc and corona state we allow 
the density and temperature of the disc and the density of 
the corona to evolve with time. We solve the equation of 
mass and angular momentum conservation: 



<9E 1 d . 1/2 d . 1/2 
at r or or 



i/E)) 



and the energy equation 
T 4 



dT dT 
dt r dr 



E 12- 10.5/3 ( ~dt +Vl ~dr' 



(18) 



(19) 



+ 



T 



PH 12 - 10.5/3 



(Q + -Q~). 



Here 



d 



Vv 



(i/Er 1/2 ) 



Er i/2 5r v'"' > ( 2 °) 

is the radial velocity in the disc while v = (2Pa)/(3pf2) is 
the kinematic viscosity. The heating term is given by Equa- 
tion [3] and the cooling term Q~ is given by Equation |H| 
while the advection is included in the energy equation via 
the radial derivatives. 

The evolution of the coronal density is given by mass 
and angular momentum conservation in the corona: 



9E cor _ 1 d 
dt r dr 



(3r ^ —— (r ^corEcor 

or 



)) + rh z 



The radial velocity in the corona is calculated as: 

cor 3 & 



cor ' 



(21) 



(22) 



E cor rV2 dr v 

with i/ cor = (2P cor a cor )/(3p cor f2) and constant viscosity pa- 
rameter in the corona a cor = 0.01. There is no need to con- 
sider the thermal evolution of the corona, since its temper- 
ature is always equal to the virial temperature and does not 
vary with time. 

We solve the above set of three time-dependent equa- 
tions using the convenient change of variables, y = 2r 1//2 
and S = yTi, at the fixed radial grid, equally spaced in y 
(see Janiuk et al. 2002 and references therein). The num- 
ber of radial zones is set to 216. After determining the 
solutions for the first 600 time steps by the fourth-order 
Runge-Kutta method, we use the Adams-Moulton predictor- 
corrector method, allowing the time-step to vary, when 
needed. 

We choose the no-torque inner boundary condition, 
£j n = Tin = for the disc. The outer boundary of the disc 
is parameterized by an external accretion rate M cx t • If this 
accretion rate is high enough, the inner disc parts gradu- 
ally heat themselves and finally end in the unstable regime, 
forcing the disc to oscillate. The boundary conditions in the 
corona are given by Equation 1151 



4 RESULTS 

4.1 Surface density and temperature evolution 

The local solutions of the accretion disc model, in the sur- 
face density vs. temperature (E- T) plane, can be calculated 
for a stationary disc in the range of accretion rates. (Alter- 
natively, on the vertical axis we can have accretion rate M 
instead of the disc temperature.) These solutions lie along 



the S-shaped stability curve, whose position on the diagram 
depends on the model parameters: black hole mass, viscosity 
and radius (c.f. Janiuk et al. 2002, Figs. 1, 3). 

Both upper and lower branches of the S-curve are vis- 
cously and thermally stable. On the lower stable branch the 
gas pressure dominates; the middle branch is unstable (ra- 
diatively cooled and radiation pressure dominant), as shown 
in detail by Pringle, Rees & Pacholczyk (1974) and Light- 
man & Eardley (1974). The upper branch is stabilized in our 
model again by the dominant gas pressure, due to the mod- 
ified viscosity law. In case of the standard viscosity (with 
constant a) this branch would be stabilized mainly by ad- 
vection, as shown in Abramowicz et al. (1988). In our case 
the advection is also taken into account, but its role is never 
dominant. 

The S-curve can also be plotted in the M — E plane, 
for any chosen disc radius. This means that the temporary 
local solutions are determined by the mean (i.e. external) 
accretion rate in the disc. Whenever the external accretion 
rate is low, so that at all the radii in the disc the local solu- 
tion sits on the lower, stable branch, the accretion proceeds 
with this rate, which is constant throughout the disc and 
constant in time. But if the accretion rate is higher than 
some critical value, M ex t > M cr it, the solutions in the in- 
nermost annuli will find themselves on the unstable branch. 
The higher M ex t, the more disc annuli will be unstable. This 
leads to the disc fluctuations, since the accretion cannot pro- 
ceed smoothly in the unstable mode. Therefore the local ac- 
cretion rate in the innermost strips changes periodically be- 
tween the lower and upper stable solutions, being no longer 
equal to A/ cxt (the accretion rate starts to depend on radius 
and time). This is displayed in the local diagrams E — T that 
are resulting form the time-dependent model. 

The exemplary stability curves of the accretion disc, 
calculated at several radii from the stationary disc model, 
are shown in Figure |3] (thin solid lines) . The thick points 
represent the subsequent solutions of the time-dependent 
model. 

The evolution of the disc on the surface density - 
temperature plane proceeds at first along the lower stable 
branch, up to the instability region. This is forced by the 
value of the external accretion rate parameter, which has to 
be large enough to drive the disc to the unstable configu- 
ration. Here we assumed M ex t = 1-5 x 10 19 g s _1 , which is 
equal to 0.45 of the Eddington rate (for black hole mass of 
10 Mq and efficiency of 1/16). The critical accretion rate in 
our model depends on whether the corona covers the disc or 
not; m cr it = 0.05 for the plain disc, while in the case of a 
disc/corona system the corona has a stabilizing role and the 
critical accretion rate is about rh cr it = 0.2. The accretion 
rates within the range 7.5 x 10 18 — 2.6 x 10 19 g s~ 1 , were 
obtained for soft states of GRS 1915+105 by Sobolewska & 
Zycki (2003). 

Firstly, in the starting, steady configuration we assumed 
the accretion rate of m ~ 1.5 x 10~ 2 of the Eddington rate 
throughout the disc. Therefore at the beginning of the subse- 
quent time-dependent calculations the model has to saturate 
at the temperatures and densities imposed by the value of 
M ox t, imposed at the outer disc radius. 

Next, the evolution proceeds in a form of loops between 
the lower and upper branches. Each loop refers to a single 
cycle of the instability, and the size of this loop depends on 
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log (£) log (S) 



Figure 3. Local evolution of the disc on the surface density - 
temperature plane, plotted for 4 values of radius: 3.53, 10.14, 
21.47 and 80.71 ftschw The thin solid line marks the stability 
curve resulting from the initial steady disc model, and the solid 
points are the time-dependent solutions. 



the location in the disc. For larger radii the loops become 
smaller and finally, in the outer disc regions (above ~ 80 — 
100J?schw), there are no instabilities and the disc remains 
stable all the time. The exact value of the maximal radius 
of the instability zone depends again on the model. In case 
of a plain disc, for the external accretion rate rh oxt = 0.45 
it is .Rmax = 100i?schw, while for m cxt = 0.56 it is i? max = 
110i?schw • In case of the disc with corona the extension of the 
unstable zone is respectively i? max = 80i?schw and 7? m ax = 
90#schw 

Note, that in the initial steady disc model we use a sim- 
plified parameterization of advection, with g a d v = TO (see 
Equation I1UI . Therefore the upper stable branch does not 
represent exactly the advective branch that results from the 
time-dependent calculations, which are based on the equa- 
tions with radial derivatives. In fact, q a dv is not constant 
throughout the disc and should depend on radius. However, 
this initial simplification does not influence our results, since 
our starting model is located on the lower, gas pressure dom- 
inated branch. Here the advection is negligible, and the sub- 
sequent time-dependent solutions match the stability curve 
perfectly. 

In Figure 2] we show the radial profiles of the surface 
density in the accretion disc, in several snapshots during 
such a loop (one instability cycle). In the minimum of a 
cycle the surface density in the inner parts of the accre- 
tion disc has a flat radial distribution. When the outburst 
starts, there appears a sharp density peak in the outer part 
of the unstable zone, which then moves outward. In the max- 
imum of the cycle this peak is accompanied by the largest 
fluctuation in the density distribution. The reason for these 
fluctuations are the viscosity and angular momentum trans- 
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Figure 4. Radial profiles of the disc surface density during the 
cycle of the evolution. The snapshots were made every 80 seconds 
and the whole cycle lasts for 625 seconds, he solid curve (1) refers 
to the minimum of the cycle (lowest disc luminosity, see also Figs. 
l6l and [SJ , the dashed curve (2) refers to the rise phase and the 
solid curve (3) refers to the maximum of the cycle. The dotted 
line (4) and other dotted lines (unnumbered) refer to the decay 
phase. 



fer changes inside this propagating "density wave". When 
the innermost radii of the disc switch to the hot state, the 
geometrical thickness of this region also rises, thus giving 
the rise to the kinematic viscosity. The increased transport 
rate results in the temporary density decay in the unstable 
zone, as the material starts to fall faster into the black hole. 
Simultaneously, at the inner edge of the disc there forms 
a temporary "bump" of material, forced by the no-torque 
inner boundary condition. 

The density fluctuation subsequently vanishes in the 
end of the cycle and during the decay phase the surface 
density in the unstable zone gradually rises, to reach the 
starting configuration. 

The corona evolves on a timescale much shorter than 
the disc. First, we investigate the corona formation in case 
of no mass exchange with the disc (Bi = B% = 0). The initial 
distribution of the surface density in the corona was flat, as 
determined by Equation 1151 When the evolution starts, the 
corona very quickly achieves its final shape - the surface 
density distribution saturates after ~125 sec. At the same 
time the disc evolves very slowly, being ready to start its 
oscillations after ~ 10 4 sec. of a gradual rise in density and 
temperature. Since there is no coupling with the corona, 
these oscillations do not influence the corona structure. 

Secondly, we proceed with the evolution in case of a 
substantial mass exchange between the disc and corona 
(B\ — B2 = 0.5). In this case, the evaporation of the disc 
accelerates the corona formation, and after ^ 100 seconds 
the coronal surface density in the maximum of the radial 
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Figure 5. Radial profiles of the corona surface density in case of 
no mass exchange with the disc (Si = B2 = 0). The profiles were 
calculated every 10 seconds and the density saturated at its final 
profile after ~ 125 seconds of the evolution. 



Figure 6. Radial profiles of the corona surface density during the 
cycle of the evolution. The snapshots were made every 80 seconds 
and the whole cycle lasts for 625 seconds. The labels are the same 
as in Fig.Q] 



distribution (around 10-Rschw) exceeds 10 g/cm 2 . Then the 
corona is further fed with material by the disc as its evolu- 
tion proceeds along the stability curve and the coronal sur- 
face density further gradually increases. The maximum sur- 
face density saturates at ~ 80 g/cm 2 after the disc reaches 
the critical point on the stability curve, which determines 
the maximum density in the disc. 

When the disc oscillations start, the corona also follows 
its time-dependent behavior. In Figure|S|we plot the surface 
density distribution in the corona during one cycle of the 
disc instability. The curves are plotted every 80 seconds in 
case of a full cycle lasting ~ 625 seconds. The solid line (1) 
is the coronal density in the minimum of the disc limit cycle, 
the dashed line (2) corresponds to the rise phase of the disc 
outburst and the solid line (3) corresponds to the maximum 
of the cycle (outburst). Immediately after the outburst the 
coronal density distribution comes back to the initial con- 
figuration (dotted lines) and remains there for the most of 
the cycle (4). Therefore during the decay phase of the disc, 
when the disc surface density gradually changes to reach a 
flat distribution, there is no change of the density in the 
corona. 

In Figure [7| we show the rate of the mass exchange be- 
tween the disc and corona, m z , during the disc instability cy- 
cle. In the minimum of the cycle the evaporation rate is very 
low and the maximal mass supply to the corona is achieved 
at r ~ 6-Rschw When the disc outburst starts, the evapora- 
tion rate grows dramatically in the middle of the unstable 
zone, while dropping to m z < at the outer edge of this 
zone. This is why the corona collapses locally at the outer 
parts, while expanding slightly in the inner parts during the 
disc outburst. In the decay phase the evaporation rate again 



becomes low, with a decaying fluctuation at the outer edge 
of the instability zone. 



4.2 Lightcurves 

The lightcurves represent the luminosity of the disc and 
corona separately. For the optically thick disc the luminosity 
is given by: 



Ldis 



4(7 

3k 



rpA 



2nrdr 



(23) 



The luminosity of the corona is calculated under the 
assumption that the corona is in the thermal equilibrium. 
The ions are heated by the viscous dissipation and either 
cool by advection or transfer their energy to the electrons, 
which in turn radiate e.g. in the Inverse Compton process. 
Therefore we have: 



Qcor 



and 
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(24) 



dr 



(25) 



In Figure|H|we show an exemplary lightcurve calculated 
for several cycles of the disc outburst. The disc limit cycle 
is very strong since M ex t = 1-5 x 10 19 g s^ 1 substantially 
exceeds the critical value. Therefore the loops marked by 
the time-dependent solutions on the E — T plane encompass 
substantial range of temperatures and densities, and the un- 
stable region of the disc has the radial extension up to about 
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Figure 7. Radial profiles of the mass exchange rate during the 
cycle of the evolution. The snapshots were made every 80 seconds 
and the whole cycle lasts for 625 seconds. The labels are the same 
as in Fig.Q] 



Figure 8. The time evolution of the disc and corona luminosity. 
Lower curve shows the luminosity outbursts of the disc, while 
the upper curve shows the simultaneous coronal fluctuations. The 
external accretion rate is Moxt = 1.5 X 10 19 g/s. 



80 -Rschw This results in regular, large amplitude outbursts 
of AlogL = 0.9 and At = 625 s. The luminosity variations 
in the corona, however of the same frequency, are not that 
strong. Also, they are sometimes anti-correlated with the lu- 
minosity changes in the disc, since at the disc rise phase the 
coronal luminosity at first drops and then rises slightly, to 
drop and rise again during the disc decay phase. Character- 
istically, between the disc outbursts the coronal luminosity 
decreases very slowly, with AlogL cor /A£ « 10 -5 , whereas 
the disc luminosity rises gradually, forming a 'wing' preced- 
ing the main outburst. 

4.3 Accretion rate 

The mass exchange rate in our prescription depends on the 
accretion rates in the disc and corona (Equation 1171 . Since 
they strongly vary with time during the cycle of the disc evo- 
lution, the mass exchange can have locally negative value. 
This has been shown in Figure Q 

In Figures l9land [TUl we show the radial profiles of the ac- 
cretion rates in the disc and corona, taken in different phases 
of the cycle: between the disc outbursts (corresponding to 
phase "1" in Fig. [7J and in the outburst peak (phase "3"). 

In the minimum of the cycle, i.e. between the outbursts, 
the accretion rate in the inner parts of the disc, where most 
of the dissipation takes place, is relatively low. It rises sub- 
stantially in the maximum of the cycle, therefore causing the 
disc luminosity outburst. The coronal accretion rate does 
not change that much, and the corresponding luminosity 
changes are not very pronounced. This is because in certain 
radii, at the outer edge of the instability strip, the radial ve- 
locity becomes negative. In consequence, the both accretion 
rates and mass exchange rate, rh z , also become negative, 



which means that some amount of material that has been 
evaporated to the corona, now goes back and sinks in the 
disc. Negative values of m z are not unphysical: stationary 
models based on the disc-corona mass exchange of Rozariska 
& Czerny (2000) predict coronal condensation at some radii 
for high accretion rates. The coronal surface density, and in 
turn the dissipation rate, is therefore reduced, even in the in- 
ner regions, as the mass is spread out over the whole corona. 
The temporary loss of the accretion rate can be estimated 
as AAf = 47rrm z Ar and it can be as high as 1.9 x 10 18 g/s. 
Therefore almost the whole rise of the coronal accretion rate 
(and luminosity), triggered by the disc outburst, is immedi- 
ately compensated by this loss due to the sinking of material 
in the strip around ~ 50-Rschw 

In addition, the accretion rate in the corona is negative 
in the outermost radii. The extension of this zone of nega- 
tive M cor depends on the outer boundary condition in the 
corona, which in our case is fixed by E cor (r ut) =l//t« 2.94. 
It implies that M COI < for r > 100i?schw, so the material 
slowly flows out from the corona at its outer radius. 

4.4 Alternative prescription for the mass 
exchange (prescription II) 

In the above used prescription for the rate of mass exchange 
between the disc and corona we expressed the total locally 
generated flux as proportional to the sum of the accretion 
rates in the disc and corona. In case of a steady-state disc 
this is equivalent to the sum of the locally generated fluxes 
by the viscous energy dissipation. However, when we con- 
sider the time evolution of an unstable disc, the radial ve- 
locity, and in turn the accretion rate, can have locally neg- 
ative values in some parts of the disc. This forces the mass 
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Figure 9. Radial profiles of the accretion rate in the disc (trian- 
gles) and in the corona (squares) between the disc outbursts. The 
time refers to the lightcurve shown in Fig. [S] 
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Figure 11. Radial profiles of the corona surface density during 
the cycle of the evolution for the case of the mass exchange rate 
given by Eq. 1261 The snapshots were made every 80 seconds and 
the whole cycle lasts for 625 seconds. The labels are the same as 
in Fig. 
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Figure 10. Radial profiles of the the accretion rate in the disc 
(triangles) and in the corona (squares) in the outburst. The time 
refers to the lightcurve shown in Fig. [S] 



prescription for the mass exchange rate can lead to different 
results of the time evolution of the disc plus corona system. 
Instead the formula for the vertical mass transfer given by 
Equation 1171 we use the following: 

= f 77T7/(r)(0.5BiQ 0O iScor— — +-B2adi 8 cJ > disc#disc) (26) 

Z GM x r ' 

This form means that the evaporation is proportional to the 
sum of the energy dissipation rates in the disc and in the 
corona, and therefore is always positive. 

In Figure tTTl we show the distribution of the coronal sur- 
face density during one cycle of the viscous-thermal insta- 
bility (cf. Fig.|SJ. Clearly, there is no local dip in the corona 
that would be caused by the negative value of m z during the 
disc outburst, as it was in the previous case. The material 
does not sink from the corona into the disc anywhere, but is 
uniformly evaporated and the rate of the evaporation is the 
highest when the disc is the most luminous. 

Similarly, the total luminosity in the corona is now cor- 
related with the disc luminosity. In Figure [T?l we show the 
lightcurves analogous to those in Figure |H] Between the out- 
bursts the coronal luminosity is slightly increasing, and the 
lightcurve exhibits regular peaks, simultaneous to the disc 
outbursts. The luminosity changes in the corona are of the 
order of AlogLcor = 0.3. 



exchange rate to locally decrease in the state of the disc out- 
burst, which is not the case for the viscously generated flux. 
In other words, the prescription for the mass exchange rate 
in the outburst of an unstable disc is no longer equivalent 
to the sum of the locally dissipated fluxes. 

Now we check whether the other, more 'conservative' 



4.5 Hard versus soft luminosity: time lags and 
luminosity-color diagram 

Since the mass exchange prescription (II) determined by 
Equation (I2(it gave more satisfactory representation of the 
time evolution with respect to the observed behaviour of 
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Figure 12. The time evolution of the disc and corona luminosi- 
ties, for the case of the mass exchange rate given by Eg, 1261 Lower 
curve shows the luminosity outbursts of the disc, while the up- 
per curve shows the simultaneous coronal outbursts. The external 
accretion rate is M ex t = 1.5 X 10 19 g/s. 



GRS 1915+105, we adopt this prescription in further anal- 
ysis. 

We calculated the cross-correlation function, as defined 
by Equation Q, for the theoretical lightcurves plotted in 
Figure 1121 Our points were separated at least by 0.039 s, 
and we found that the corona lightcurve lags the disc by 
1.16 seconds. This is shown in Figure ESI The time lag m 
this case is equal to about 0.5% of the duration of an out- 
burst. We note that the outburst duration depends on the 
external accretion rate M ex t and the lag depends on the 
coronal viscosity parameter a, so we could obtain a range of 
values for different model parameters. 

In Figure ITU we show the luminosity-color diagram for 
one of the observations of GRS 1915+105, that shows the 
pronounced outbursts (class p, cf. Tabled - On the horizon- 
tal axis we plot the soft X-ray flux in the range 1.5-6 keV, 
and on the vertical axis we plot the X-ray color, i.e. the ra- 
tio of the hard to soft fluxes, -F , 6.4-l4.6fceV r /-Fi.5-6fceV ■ During 
its evolution, the source follows a characteristic track, in the 
form of a loop on this diagram. In the lower region (the "ba- 
nana" shape) the soft luminosity is low, and the X-ray color 
is soft. The upper region ("island") shape is characterized by 
substantial luminosity both in hard and soft X-ray bands. 

In Figure H31 we show the result of our modeling. On the 
horizontal axis of the theoretical luminosity-color diagram 
we plot the disc luminosity, while in the vertical axis we plot 
the ratio of the corona to disc luminosities. The "banana" 
shape is reproduced in our model quite well. However, we 
do not obtain the other distinct region on this diagram (the 
"island"). In our simulations the source moves only along 
the track shown in Fig. 1151 and after reaching the end of the 
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Figure 13. The normalized cross-correlation function of the disc 
and corona lightcurves. The main peak is shifted to -1.16 s, which 
means that the corona lightcurve lags the disc one by 1.16 s. 

"banana" region immediately jumps to the top-left corner 
of the diagram. 



5 DISCUSSION 

The main simplification of our model was the assumption 
that the temperature in the corona is constant and equal 
to the virial one. Therefore the corona is heated only via 
increase of its surface density and is not radiatively coupled 
to the disc. We do not consider here any specific radiative 
processes that are important in cooling of the coronal gas, 
e.g. Comptonization. 

On the other hand, this approach lets us solve the set of 
three time-dependent equations for the disc plus corona evo- 
lution, accompanied by the formula for the mass exchange 
rate, and avoid numerical problems with the fourth equa- 
tion for the corona temperature. Using a finite time-step we 
are able to follow the time evolution of the system as long 
as it saturates in a quasi-stationary state, with oscillations 
of constant amplitude and duration. We qualitatively and 
quantitatively check, how the behaviour of the underlying, 
unstable disc affects the corona and vice versa. 

The full, two- or three dimensional treatment to 
the global long-term evolution of an accretion disc (or 
disc/corona system) is very complex. Up to now, the so- 
phisticated 2-D and 3-D accretion disc simulations (e.g. Agol 
et al. 2001; Turner et al. 2003) either treated the problem 
locally, or did not include radiative cooling. Turner (2004) 
made the first attempt to compute the 3-D simulation in 
the flux-limited approximation; in this case however the 
plasma temperatures are still much too low to form the hot 
corona. Since the simulations are not able to reach the vis- 
cous timescale of the disc, the long-term evolution of an 
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Figure 14. The luminosity-color diagram of the two observed 
lightcurves, shown in Fig. On the horizontal axis is the soft 
X-ray flux in the range 1.5-6 keV, and in the vertical axis is the 
ratio of the hard to soft fluxes, .F6.4-l4.6fcev/i ; l.5-6fceV- 
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Figure 15. The luminosity-color diagram of the two theoretical 
lightcurves, shown in Fig. 1121 On the horizontal axis is the soft 
X-ray (disc) luminosity, and in the vertical axis is the ratio of the 
hard to soft (corona to disc) luminosities. 



X-ray source cannot be followed and any observational con- 
sequences of such a model are difficult to be checked. There- 
fore the simple approach presented in this paper can still be 
valuable. 



5.1 Viscosity parameterization 

In our choice of the viscosity law we follow the paper of 
Nayakshin, Rappaport & Melia (2000). The motivation of 
these authors was mostly observational: they needed a vis- 
cosity law which allows for a disc instability at the inter- 
mediate accretion rates but provides the stable solution at 
low and high accretion rates. Moreover, the stable upper 
branch should appear at accretion rates smaller or compa- 
rable to the Eddington rate. The standard aPtot viscosity 
law of Shakura & Sunyaev (1973) satisfies in a natural way 
all requirements but the last one. Upper stable branch in 
such solution forms due to advection (i.e. the "slim disc"; 
Abramowicz et al. 1988), and the efficient advection develops 
only for accretion rates much larger than the Eddington val- 
ues, definitely too high too account for the time behaviour 
e.g. of the microquasar GRS 1915+105. What is more, in 
most cases the unstable branch should cease to exist at all 
since most of the X-ray binaries accreting at high (but sub- 
Eddington) rates are quite stable and well described by the 
classical disc (Gierliriski & Done 2004). 

Therefore, observations tell us that some modifications 
are absolutely needed. Either we must postulate very strong 
outflow (which effectively cools the disc), or a modification 
of the dissipation law itself. 

We have much better understanding and numerous ob- 
servational constraints for the viscosity law when the gas 
pressure dominates. 3-D MHD simulations of the magneto- 
rotational instability (MRI) well explain the nature of the 
angular momentum transfer and the rough value of the vis- 
cosity parameter. Still, some ad hoc modifications, like ad- 
ditional dependence on the disc thickness, are used to model 
the development of the ionization instability responsible for 
outbursts in cataclysmic variables and numerous X-ray tran- 
sients (e.g. Cannizzo et al. 1995). 

The theoretical background of the viscosity law is still 
quite poor in case when radiation pressure is important. Sev- 
eral authors in the past suggested that actually the whole 
idea of aPtot scaling is inappropriate, and instead the aP gas 
should perhaps be used even if the radiation pressure dom- 
inates (Lightman & Eardley 1974; Stella & Rosner 1984). 
Sakimoto & Coroniti (1989) argued that the magnetic field 
will be expelled from the radiation pressure dominated disc 
by the buoyancy so the effectiveness of the angular mo- 
mentum transfer must decrease with an increase of radia- 
tion pressure. Other authors argued that whenever radiation 
pressure dominates, new kinds of instabilities develop which 
may modify the disc structure considerably. Gammie (1998) 
discussed the photon bubble instability, and Ruszkowski & 
Begelman (2003) argued this instability leads to the disc 
dumpiness which in turn influences the radiative transfer. 

We expect some improvement of the viscosity model- 
ing with the future development of MHD simulations of ra- 
diation pressure dominated discs. So far, two such simula- 
tions were performed. 2-D computations of Agol et al. (2001) 
lasted only for a fraction of the disc thermal timescale and as 
a result the disc collapsed to gas pressure dominated state. 
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A 3-D simulation of the disc which reached the thermal 
balance and lasted for about 8 thermal timescales was re- 
cently completed by Turner (2004). The disc did not achieve 
the full stability and showed long-lasting variations by a 
factor of a few but the instability was not as violent as pre- 
dicted by the standard aP to t mechanism (Szuszkiewicz & 
Miller 1998; Janiuk et al. 2002). The magnetic field lines 
were indeed partially expelled from the disc interior. In this 
simulation, the initial Shakura & Sunyaev (1973) state with 
a = 0.01 evolved to a complex state with mean dissipation 
level equivalent to a = 0.0013. 

Therefore, the modified viscosity law, given by Eq. 
seems to be a plausible option in case of the accretion disc. 
Since in the corona the radiation pressure does not con- 
tribute to the total pressure, the classical viscosity parame- 
terization of the corona with constant a is also justified. 

5.2 Comparison with observations 

Our time dependent model of the non-stationary accretion 
onto a black hole gives at least a partial explanation of the 
complex variability of the microquasar GRS 1915+105. Di- 
rect evidence for variable accretion rate in this microquasar 
comes from the spectral analysis (Migliari & Belloni 2003), 
and a plausible mechanism that leads to the local accre- 
tion rate variations in the inner parts of a disc, giving the 
outbursts of appropriate amplitudes and durations, is the 
radiation pressure instability. 

The microquasar is unique because of its large ampli- 
tude, regular outbursts. Such a behaviour is not observed 
in other Galactic X-ray sources, probably because the Ed- 
dington ratio in this microquasar is higher than in other 
black hole systems (Done, Wardzinski & Gierlihski 2004). 
Our model, at least qualitatively, explains this phenomenon. 
The radiation pressure instability leads to the strong out- 
bursts only if the accretion rate is higher than some critical 
value. Also, the corona has a stabilizing role and can sup- 
press the disc oscillations completely or make them less pro- 
nounced. Quantitatively, however, the reliable determina- 
tion of the critical accretion rate above which the outbursts 
occur, would require a more detailed modeling of viscosity 
within the disc and essentially the 2-dimensional calcula- 
tions. This is beyond the scope of the present work. 

The physical coupling between the disc and corona via 
the mass exchange leads in a natural way to the regular time 
delays between the disc (soft X-ray) and corona (hard X- 
ray) emission. The lag of the order of 1 s is comparable with 
the viscous timescale in the corona and is required by this 
hot flow to adjust to the variable conditions in the underly- 
ing disc. Such lags are present in some observed lightcurves 
of the microquasar, i.e. these exhibiting outbursts that are 
possibly connected with the radiation pressure instability. 

The strong variability of GRS 1915+105 manifests it- 
self also by characteristic tracks on the color-color and 
luminosity-color diagram. The bottom-left corner of the 
luminosity-color plot is occupied by the source during its 
gradual rise (the "wing"), preceding the outburst state. 
This, so called here "banana" shape is well modeled in our 
calculations. While the disc luminosity (soft X-ray flux) is 
gradually rising before the outburst, the rise in the coronal 
luminosity is much flatter, leading therefore to the decreas- 
ing hard-to-soft flux ratio (color). On the other hand, in the 



outburst peak, when the soft luminosity is the highest, the 
color starts rising. Just after the outburst there appears a 
somewhat flat maximum in both hard and soft X-rays, lead- 
ing to the other distinct region occupied by the source on 
the luminosity color-diagram: so called here "island" shape. 
This is not present in our model calculations, since our out- 
burst maximum is a sharp peak. Therefore we do not obtain 
a state of both high disc and corona fluxes, which in conse- 
quence would give the high color in the disc emission peak. 
Instead, the corona lags the disc emission and its maximum 
corresponds to the already decaying state of the disc. 



6 CONCLUSIONS 

We presented the first results of the time-dependent calcu- 
lations of thermal- viscous evolution of accretion disc that is 
coupled to the corona by the mass exchange. Contrary to 
the accretion disc, the corona is stable, since there is only 
gas pressure included in its equation of state. However, it 
is also to a certain extent subject to similar changes as the 
accretion disc during its instability cycle, since it actively 
responds to the behaviour of the underlying disc. 

The main conclusions of the present work are as follows: 

• The luminosity outbursts in the disc are correlated with 
similar outbursts in the corona if the mass exchange rate is 
proportional to the sum of the locally dissipated fluxes ( 
prescription II). The assumption that it is proportional to 
the sum of the local accretion rates (prescription I) is equiv- 
alent, but only during the disc quiescence. In the outburst 
this prescription leads to the luminosity dips in the corona 
rather than outbursts. Specifically, the lightcurves of the disc 
and corona are anti-correlated, with two dips in the corona 
corresponding to the peak of the disc luminosity. 

• The luminosity profile between the outbursts in the 
corona is gradually decreasing, while in the disc it is grad- 
ually rising (prescription I). Alternatively, the prescription 
II leads to the slow rise of the coronal flux, correlated with 
the steep rise of the disc. 

• The coronal outbursts are of much smaller amplitude 
than the disc ones, but of similar duration. 

• The outbursts in the corona lag the disc evolution by 
~ 1 second, which is in good agreement with observations 
of the microquasar GRS 1915+105. 

• The outer boundary condition in the corona determines 
if the material flows out at the outer edge. For small surface 
density E out ~ 3 g/cm 2 the accretion rate in the corona is 
negative above ~ lOCRschw 

• If the mass exchange between the disc and corona is 
proportional to the sum of the local accretion rates (pre- 
scription I), the material may also locally sink into the disc 
in the region around 50i?sohw 

Our model is only in part able to reproduce the tracks 
of the system on the luminosity-color diagram, having a "ba- 
nana" shape without the "island" one. This should be mod- 
eled in more detail, including also the spectral analysis. In 
addition, the tracks on the color-color diagram, not analyzed 
here, in some observed sources may have a form of hysteresis 
(Maccarone & Coppi 3002; Zdziarski et al. 2004). Since the 
phenomenon is connected with much longer timescales, than 
those considered here, the underlying instability mechanism 
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should be different. A plausible one seems to be the partial 
hydrogen ionization instability, which operates essentially in 
the same way, causing the disc to oscillate between the cold 
and hot phases. This is the subject of our future investiga- 
tions. 
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